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Interpolating  Functions  on  Lines  in  3-Space 


Martin  Peternell  and  Helmut  Pottmann 


Abstract.  Given  straight  lines  Li,  i = 1, . . . , N,  in  Euclidean  3-space 
with  associated  function  values  /,,  we  study  the  interpolation  problem  of 
constructing  a smooth  real  valued  function  F which  interpolates  values  fi 
at  given  data  lines  L The  function  F shall  be  defined  on  the  entire  set 
of  lines  or  at  least  on  lines  contained  in  a domain  of  interest  in  3-space. 


§1.  Introduction 

The  problem  of  constructing  an  interpolating  function  F for  data  lines  L,  and 
corresponding  function  values  /,  is  a scattered  data  interpolation  problem  in 
the  set  of  lines  C in  Euclidean  3-space  E3 . 

A variety  of  solutions  of  scattered  data  interpolation  problems  for  data 
points  Xi  € U with  U = !Rn  or  U C IR"  are  known,  see  [3].  Extensions  to 
spheres  and  other  surfaces  in  1R3  are  described  in  [2]  and  references  therein. 

Scattered  data  interpolation  on  lines  is  quite  different,  since  the  set  of  lines 
C is  not  a Euclidean  space.  It  is  a result  of  classical  geometry  that  the  set  of 
lines  C of  projective  extension  P 3 of  Euclidean  3-space  E3  is  a 4-dimensional 
quadratic  variety  'n  projective  P5.  Thus,  the  general  formulation  of  the 
problem  is  as  follows:  Construct  a function  F : M4  — ► IR  interpolating  val- 
ues fi  to  corresponding  data  lines  L,.  For  practical  purposes  it  is  sufficient 
to  construct  (or  represent)  functions  on  subsets  of  M4  which  correspond  to 
domains  of  interest  in  E3,  containing  all  data  lines. 

The  solution  presented  here  will  be  the  following.  We  restrict  to  specific 
four-dimensional  subsets  Co  of  . These  subsets  possess  parametrizations 
IR4  — > Co  with  the  property  that  distances  between  lines  in  Co  are  induced 
by  special  positive  quadratic  forms  in  IR4.  This  fact  allows  us  to  apply  well- 
known  methods  in  IR4  to  solve  interpolation  (or  also  approximation)  problems. 

Applications  include  light  field  rendering  in  computer  graphics  [4].  Con- 
sidering motion  planning  in  robotics,  the  method  applies  to  represent  a dis- 
tance function  of  robot  arms  (lines)  to  obstacles.  The  first  motivation  for 
studying  functions  on  lines  came  from  five  axis  milling.  There,  the  question 
occurs  of  how  to  represent  axis  positions  (lines)  of  the  cutting  tool. 
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§2.  Lines  in  Space 

An  oriented  line  L in  Euclidean  3-space  E3  is  determined  by  a point  p and  a 
unit  direction  vector  1 (||1||  = 1).  Together  with  the  moment  vector 

i = pxi,  (i) 

we  obtain  a representation  of  L by  a sixtuple 

L = (1;I)  = {h,h,h\Uih,h)-  (2) 

These  Ij’s  are  called  normalized  Pliicker  coordinates  of  L.  By  (1),  these  coor- 
dinates are  not  independent,  but  satisfy  the  Pliicker  relation 


1 • 1 — l\li  + ^2^5  + ^3^6  — o.  (3) 

Substituting  1 by  —1  leads  to  coordinate  vector  — L which  defines  the  same  line 
but  with  opposite  orientation.  To  get  more  information  about  the  structure 
of  lines  in  space,  it  is  necessary  to  study  the  set  of  lines  C in  the  projective 
extension  P3  of  E3. 

E3  is  extended  to  P3  by  adding  points  and  lines  at  infinity.  Using  the 
analytical  model  IR4,  points  in  P3  are  one  dimensional  subspaces  of  R4.  Thus, 
we  will  use  the  following  notation  for  points  in  P3, 

(xo,xi,X2,  £3)R  :=  (Ax0,  • • • , Ax3),  A £ R. 

Let  w : xo  = 0 be  the  plane  at  infinity.  We  write  briefly  (xo,  x)R,  with  x £ R3 
for  points  in  P3.  The  transition  from  homogeneous  to  Cartesian  coordinates 
is  given  by 

/ \ / *^1  ^2  \ 

32,33)  ( — , — , — 

30  3q  Xo 

which  is  obviously  only  possible  for  points  not  at  infinity. 

A line  L in  P3  usually  is  spanned  by  two  points  (po,p)R  and  (go,q)R- 
Homogeneous  Pliicker  coordinates  are  obtained  by 

L = (ii,...,Z6)  = (poq-9oP,pxq).  (4) 

If  we  substitute  (po,p)  by  A(po,p),  we  get  AL  such  that  the  /,’s  are  only 
determined  up  to  a scalar  multiple.  This  proves  homogeneity  of  L. 

If  L is  not  in  w,  the  relation  to  definition  (2)  is  obtained  as  follows. 
Let  (p0,p)R  be  a proper  point  on  L such  that  we  can  switch  to  Cartesian 
coordinates  p by  letting  po  = 1.  Further,  let  (q0,  q)R  be  the  intersection 
point  ufl  L which  implies  qo  = 0.  Inserting  this  into  (4)  gives  (2)  up  to  a 
normalization  of  the  direction  vector  q = 1 of  the  line  L. 

If  L is  in  w,  its  Pliicker  coordinates  are  (o,  a)R  with  o = (0, 0, 0)  and 
some  not  vanishing  vector  a.  We  can  interprete  L as  the  line  at  infinity  of  a 
pencil  of  parallel  planes  a • x = c,  with  c £ R.  All  these  planes  possess  a as 
normal  vector. 
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Fig.  1.  Stereographic  projection  of  a hyperboloid  Q. 

Since  L and  AL  define  the  same  line  in  P3,  homogeneous  Pliicker  coordi- 
nates (4)  define  points  LR  in  P5.  But  only  those  6-tuples  (x\, . . . , 26)R  are 
Pliicker  coordinates  of  a line  X in  P3,  which  satisfy 

X1X4  -(“  X2X5  -j-  x$xq  — 0. 

This  quadratic  variety  is  called  the  Klein  quadric  , where  upper  and  lower 
indices  denote  dimension  and  degree  of  this  variety.  The  maximal  dimension 
of  its  subspaces  is  2.  It  is  a point  model  of  the  set  of  lines  £ of  P3.  The 
bijection 

7 : £ ->  M| 

from  lines  L C P3  to  points  LR  of  M2  is  called  Klein  mapping. 

The  image  points  (o,  a)R  of  lines  at  infinity  lie  in  the  plane  Eu  : x\  = 
x2  — X3  = 0 which  is  entirely  contained  in  M2 . All  lines  passing  through 
the  origin  O = (1,0, 0, 0)R  have  Pliicker  coordinates  L = (1,  o).  This  can  be 
checked  by  letting  p = (0, 0, 0)  in  formula  (1).  The  corresponding  image  points 
in  P5  lie  in  the  plane  E0  : X4  = X5  = xq  = 0.  In  general,  all  lines  through  an 
arbitrary  point  in  P3  possess  7-images  which  lie  in  a 2-dimensional  subspace 
of  M2  ■ The  same  holds  for  lines  contained  in  an  arbitrary  plane  in  P3.  Thus, 
M2  contains  two  3-parametric  families  of  2-dimensional  subspaces. 

We  emphasize  that  C and  £ are  not  Euclidean,  affine  or  projective  spaces. 

Local  coordinates  of  lines 

We  have  seen  that  C is  isomorphic  to  M2  — Ew,  where  Eu  consists  of 
image  points  of  all  lines  at  infinity.  Let  T be  the  tangent  hyperplane  of  M| 
at  a point  Z and  let  r = M|  fl  T.  It  is  known  that  r is  the  7-image  of  all  lines 
intersecting  the  line  L = Z 7-1. 

Lemma  1.  M2  — r = A4  is  an  affine  space. 

Proof:  This  lemma  is  a result  of  classical  geometry,  and  is  proved  by  stere- 
ographic projection.  Let  Q be  a regular  quadric  in  P".  Let  Z be  a point  in 
Q and  T its  tangent  hyperplane,  see  Figure  1.  Further,  consider  E to  be 
a hyperplane  in  Pn , not  incident  with  Z.  The  intersection  r = Q fl  T is  a 
quadratic  cone  with  vertex  Z.  The  intersection  e = E fl  T is  a hyperplane 
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Fig.  2.  Local  coordinates,  distance  function. 

in  E.  This  says  that  E — e is  an  affine  space.  The  stereographic  projection 
cr  : Q — t —>  E — e with  center  Z is  bijective  and  maps  points  P £ Q — t to 
points  P'  in  affine  space  E — e.  □ 

Figure  1 shows  a low  dimensional  example.  Q is  a hyperboloid,  and  r is 
a pair  of  lines.  Planes  E and  T are  parallel  such  that  e is  at  infinity. 

We  come  back  to  line  geometry  and  the  Klein  quadric  M\-  Let  Z = 
(0, 0, 0, 0, 0, 1)1R  be  the  center  of  a stereographic  projection.  It  is  the  7-image 
of  the  line  at  infinity  which  is  determined  by  horizontal  planes  2 =const.  with 
normal  vector  (0,0,1).  The  tangent  hyperplane  T at  Z with  respect  to 
is  given  by  the  equation  £3  = 0.  The  exceptional  set  r = M.j  D T consists 
of  7-images  (ly,  I2, 0; . . .)1R  of  all  horizontal  lines.  Lemma  1 says  that  all  non- 
horizontal lines  form  an  affine  space  A4. 

Consider  two  horizontal  planes  Eo  : z = 0 and  E\  : z = 1.  The  inter- 
section points  go  = (<?i,g2,0)  and  gi  = (93,94,1)  of  a line  G and  planes  Ey 
(Figure  2)  define  a parametrization  of  all  non-horizontal  lines  by 

1R4  = IR2  x 1R2  — > £ 

, s r (5) 

(91,92,93,94)  ^ G. 

P liicker  coodinates  of  G are  G = (93-91,94  -92, 1;  92,  ~9i , £7i94  - 9293)-  The 
stereographic  projection  with  center  Z onto  £6  = 0 gives 

G'  = (93  - 9i>94  - 92, 1;92,  — 9i > 0) . 

This  equals  (5)  up  to  a linear  mapping.  Hence,  the  mapping  (5)  from  non- 
horizontal lines  to  points  in  IR4  is  geometrically  equivalent  to  a stereographic 
projection  of  M4  — r. 

Distance  function  of  lines 

For  practical  purposes,  it  is  sufficient  to  consider  distances  of  lines  within  a 
domain  of  interest.  To  specify  this  domain,  we  will  consider  only  lines  which 
enclose  an  angle  < cf> 0 with  a fixed  unit  vector  z.  The  unit  direction  vector  g 
of  such  a line  G satisfies 


g • z > cos  rf>  0 ■ 
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We  have  chosen  a Cartesian  coordinate  system  with  z as  third  axis.  Further, 
we  will  consider  only  segments  of  lines  between  two  planes  Eq,Ei,  bounding 
the  domain  of  interest.  This  is  motivated  by  the  fact  that  we  are  interested  in 
particular  in  distances  between  points  lying  between  those  planes.  Let  gj,hj 
be  intersection  points  of  lines  G,H  with  Ei,  and  consider  points  x,  y on  G 
and  H,  respectively, 

x(t)  = (1 -t)g0 +tgi, 
y(t)  = (1  — £)h0  + thi. 

The  square  of  a useful  distance  between  lines  G,  H within  the  above  domain 
of  interest  is  defined  by 

d(G,  H)2  = ||x(t)  - y(t)\\2dt  ^ 

= (go  - h0)2  + (gi  - hi)2  + (g0  - h0)  • (gi  - hi). 

It  measures  horizontal  distances  between  corresponding  points  x,  y of  G,  H . 
We  will  not  distinguish  between  a line  X and  its  coordinate  vector  X — 
(aq,  X2,  X3,  £4)  in  1R4  according  to  parametrization  (5).  Formula  (7)  is  a posi- 
tive definite  quadratic  form  in  H4  with  the  following  coordinate  representation 

(X,  X)  = x\  + x2  + £3  + x\  + 2:1X3  + X2X4. 

Remark  2.  These  distances  differ  from  orthogonal  distances  (from  a point  to 
a line  G)  only  by  a factor  < cos  4>o-  So,  taking  (j> 0 relatively  small  will  control 
the  difference  between  these  distances  and  the  Euclidean  distances  in  E3. 

Summary  3.  The  restriction  to  specific  subsets  Co  of  line  space  allows  para- 
metrizations  ]R4  — > £q.  A positive  definite  quadratic  form  in  1R4  serves  to 
define  distances  between  lines  in  a useful  manner. 

Choice  of  local  coordinates 

Distance  d is  not  invariant  under  motions  in  E3,  but  depends  on  the  choice 
of  z and  planes  Eq,Ei.  Consider  oriented  lines  Li,i  = 1 with  unit 

direction  vectors  1;.  Assume  that  lj  • 1 k < C.  This  expresses  that  the  angle 
between  any  two  lines  is  bounded  by  arccos(C).  A good  choice  for  the  vector 
z can  be  computed  as  solution  of  a regression  problem.  Assuming  ||lj||  = 1, 
we  want  to  maximize 

f>-3)2  (8) 

i= 1 

over  all  unit  vectors  z.  Maximizing  the  quadratic  form  (8)  under  the  quadratic 
side  condition  z • z = 1 leads  to  an  eigenvalue  problem  in  1R3.  Thus,  we 
found  a possibility  to  construct  z with  respect  to  a set  of  lines  Li.  Planes 
Eo,  Ei  perpendicular  to  z bounding  the  domain  of  interest  have  to  be  chosen 
depending  on  the  problem.  In  this  sense  we  can  say  that  the  coordinate  system 
is  connected  with  the  problem  in  an  invariant  way. 
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Fig.  3.  Definition  of  domains. 

If  direction  vectors  1*  of  lines  L,  are  distributed  over  a whole  hemisphere  or 
more,  we  have  to  split  the  set  of  lines  into  subsets  and  perform  the  construction 
of  coordinate  systems  for  the  subsets.  Remark  2 gives  information  about  the 
deviation  of  distances  compared  to  usual  distances  in  E 3. 

§3.  Representation  of  Functions  on  £ 

Given  N lines  L,  with  corresponding  function  values  /,,  we  would  like  to 
compute  a function  F : £ — » IR  with  F(Li)  = fi . This  is  a scattered  data 
interpolation  problem  on  £ (or  )■  With  help  of  local  parametrizations  we 
obtain  scattered  data  interpolation  problems  on  1R4.  The  given  algorithm 
consists  of  three  steps. 

1)  Find  a covering  {Uj,j  = of  £ with  domains  Uj  which  are 

parametrized  over  JR4.  Decide  the  membership  of  lines  and  domains. 

2)  Compute  partial  solutions  Fj  of  the  interpolation  problem  for  all  domains 

Uj- 

3)  Merge  all  partial  solutions  Fj  in  a global  solution  F with  required  conti- 
nuity. 

First  of  all  we  want  to  find  a covering  of  lines  Li  by  domains  Uj  with  1 < 
j < M . We  choose  M unit  vectors  z j and  real  numbers  Rj  which  serve  as 
centers  and  spherical  radii  of  caps  of  the  unit  sphere  S2.  These  caps  determine 
domains  Uj  in  the  following  way.  A line  L belongs  to  Uj  if  and  only  if 

1 ■ Zj  > cos  Rj 

holds  for  its  direction  vector,  see  Figure  3.  Clearly,  L can  be  contained  in  more 
than  one  domain.  We  determine  the  membership  of  all  lines  L;  for  domains 

Uj- 

In  a second  step  we  compute  partial  solutions  Fj  of  the  interpolation 
problem  for  each  domain  Uj . This  is  done  by  letting 

N, 

Fj(X)  = Y,«ikBk{X), 

k= 1 
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where  Nj  shall  be  the  number  of  lines  Li  belonging  to  domain  Uj.  X = 
(xi,X2,xs,Xi)  € 1R4  is  coordinate  vector  of  a line  X according  to  parametriza- 
tion  (5).  Bk(X)  are  (for  instance)  radial  basis  functions  and  depend  only  on 
the  distance  d(X,  Lk )•  The  coefficients  ajk  are  solutions  of  linear  systems.  The 
problem  of  regularity  of  such  systems  dependent  on  the  type  of  basis  function 
is  solved  in  [5].  So  we  get  partial  solutions  Fj  valid  in  domains  Uj. 

In  the  last  step  we  have  to  merge  all  partial  solutions  to  a unique  one. 
This  can  be  done  by  forming  a weighted  sum 

M 

F(X)  = Y,wi(X)Fi(X)- 

j= i 

The  weights  can  be  chosen  as 

= (1  ~ arccos(x  • mj)/Rj)r+ 

3 Ei=i(1-arccos(x'm0/-Ri)+’ 

where  m_,  and  Rj  are  center  and  radius  of  the  spherical  cap  which  defines  Uj 
and  x denotes  the  normalized  direction  vector  of  the  line  X.  The  notation 
(q)T+  expresses  that  Wj(X)  is  positive  in  the  interior  of  Uj  and  is  zero  outside. 
This  says  that  (<?)+  = qT  for  positive  q,  and  (q)T+  = 0 otherwise. 

Weights  wj(X)  are  in  the  differentiability  class  Cr_1.  If  partial  solutions 
Fj  possess  the  same  smoothness,  then  also  F is  in  Cr~1. 

§4.  Visualization  of  Functions  on  Lines 

Since  the  dimension  of  C is  four,  visualization  of  function  values  is  an  advanced 
topic.  In  general,  displaying  functions  on  low  dimensional  subsets  seems  to  be 
promising.  We  decided  to  choose  several  bundles  of  lines  for  evaluation  and 
want  to  describe  two  methods  of  visualization. 

We  choose  an  appropriate  number  of  points  v , within  the  domain  of 
interest,  and  evaluate  F at  sufficiently  many  lines  passing  through  vertices 
Vj.  Let  Fmax  be  an  (existing!)  upper  bound  of  the  absolute  function  values. 
Consider  lines  Lij  with  function  values  F(Lij)  = fij  passing  through  vertex 
Vj.  Assume  that  Ly  are  oriented  lines.  Displaying  the  star-shaped  surfaces 

Pi  = Vi  + (1  + 

for  all  chosen  vertices  v*  is  one  possibility  to  visualize  function  values.  If 
function  values  for  L and  — L are  equal,  the  p,  will  be  centrally  symmetric 
surfaces.  For  functions  on  nonoriented  lines,  we  will  use  both  direction  vectors 
ly  and  — ly  for  the  definition  of  py  and  assign  the  same  function  value  /y  to 
them.  Thus  we  always  get  centrally  symmetric  surfaces.  Figure  4 shows  an 
interpolant.  The  test  function  is  a function  of  the  distances  between  lines  Lj 
and  points  (not  displayed). 

For  the  second  method  we  use  spheres  Si,  centered  at  vertices  V;.  All 
lines  L^  of  the  bundle  v,  with  constant  function  values  form  a cone  C with 
vertex  v,.  Intersecting  these  cones  C(ct)  for  several  constants  c*  gives  level 
curves  on  spheres  Si  (not  displayed). 
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Fig.  4.  Visualization  of  functions  on  lines. 
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